\section{Finite Element Formulation of Elastic Solid Dynamics}

Here we gave a brief derivation of using finite element method to discretize elastic solid dynamic governing equation. In short, the motion of elastic solid can be described by the following equations,

\begin{equation}\nonumber
\begin{cases}
\rho \ddot{\textbf{u}}(\textbf{x}) + \nabla\cdot\boldsymbol{\sigma}(\textbf{x}) + \textbf{f}_{b}(\textbf{x}) = \textbf{0}, \textbf{x} \in \Omega\\
\textbf{u}(\textbf{x}) = \textbf{u}_d(\textbf{x}), \textbf{x} \in \partial\Omega_D\\
\boldsymbol{\sigma}(\textbf{x})\textbf{n}(\textbf{x}) = \textbf{f}_{s}(\textbf{x}), \textbf{x} \in \partial\Omega_N
\end{cases}
\end{equation}

Here \textbf{x} represents position in material space $\Omega$, \textbf{u} is a displacement field which can be taken as a mapping from material space to world space, $\boldsymbol{\sigma}$ is stress tensor, $\textbf{f}_b$ is body force while $\textbf{f}_s$ is surface traction, \textbf{n} is normal field or Gauss map. Finally, $\partial\Omega_D$ and $\partial\Omega_N$ represents Dirichlet boundary and Neumann boundary respectively. 

This is the dynamic governing equation. To make our life easier, we omit the damping term and incompressibility condition which can be turned into a constrained dynamics problem and solved by constructing a KKT matrix with pressure as Lagrange multiplier.

Another simplification is we assume linearity in both geometry and material, which can be phrased as following

\begin{equation}\nonumber
\begin{cases}
\boldsymbol{\epsilon} = \textbf{L}(\textbf{u})\\
\boldsymbol{\sigma} = \mathcal{C} \colon \boldsymbol{\epsilon}
\end{cases}
\end{equation}

Here $\textbf{L}(\cdot)$ is a linear first order differential mapping and $\mathcal{C}$ is a forth order elasticity tensor.

Now that we have our problem clearly defined, we can start discretization using FEM. We assume that both our trial space and test space are $H^{1}(\Omega)$. Then we apply Galerkin projection to the dynamical equation and obtain its weak formulation

\begin{equation}\nonumber
\begin{split}
&\int_{\Omega} <\textbf{v}(\textbf{x}), \rho \ddot{\textbf{u}}(\textbf{x}) + \nabla\cdot\boldsymbol{\sigma}(\textbf{x}) + \textbf{f}_{b}(\textbf{x})> d\Omega = 0, \forall \textbf{v}\\
\Rightarrow &\int_{\Omega}<\textbf{v}, \rho \ddot{\textbf{u}}>d\Omega + \int_{\Omega} <\textbf{v}, \nabla\cdot\boldsymbol{\sigma}>d\Omega + \int_{\Omega} <\textbf{v}, \textbf{f}_{b}>d\Omega = 0\\
\Rightarrow &\int_{\Omega}<\textbf{v}, \rho \ddot{\textbf{u}}>d\Omega + \int_{\partial\Omega} <\textbf{v}, \boldsymbol{\sigma}\textbf{n}>dA - \int_{\Omega}\nabla\textbf{v} \colon \boldsymbol{\sigma}d\Omega + \int_{\Omega} <\textbf{v}, \textbf{f}_{b}>d\Omega = 0\\
\Rightarrow  &\int_{\Omega}<\textbf{v}, \rho \ddot{\textbf{u}}>d\Omega + \int_{\partial\Omega_N} <\textbf{v}, \boldsymbol{\sigma}\textbf{n}>dA - \int_{\Omega}\nabla\textbf{v} \colon \boldsymbol{\sigma}d\Omega + \int_{\Omega} <\textbf{v}, \textbf{f}_{b}>d\Omega = 0\\
\Rightarrow &\int_{\Omega}<\textbf{v}, \rho \ddot{\textbf{u}}>d\Omega + \int_{\partial\Omega_N} <\textbf{v}, \textbf{f}_s>dA - \int_{\Omega}\nabla\textbf{v} \colon \boldsymbol{\sigma}d\Omega + \int_{\Omega} <\textbf{v}, \textbf{f}_{b}>d\Omega = 0\\
\Rightarrow &-\int_{\Omega}<\textbf{v}, \rho \ddot{\textbf{u}}>d\Omega + \int_{\Omega}\nabla\textbf{v} \colon \boldsymbol{\sigma}d\Omega = \int_{\partial\Omega_N} <\textbf{v}, \textbf{f}_s>dA + \int_{\Omega} <\textbf{v}, \textbf{f}_{b}>d\Omega
\end{split}
\end{equation}

Now we move from infinite dimension to finite dimension by introducing subspace $H^{1}_{h}(\Omega)$ of $H^{1}(\Omega)$, which can be spanned by a finite number of basis functions $\phi_i$ selected from $H^{1}(\Omega)$. Thus the true solution can be approximated by the following

\begin{equation}\nonumber
\textbf{u} = \sum_{i} u_i \boldsymbol{\phi}_i
\end{equation}

By plugging the above expansion into the weak formulation, we can have

\begin{equation}\nonumber
\begin{split}
&-\int_{\Omega} <\textbf{v}, \rho\sum_{i}\ddot{u}_i\boldsymbol{\phi}_i>d\Omega + \int_{\Omega}\nabla\textbf{v} \colon \boldsymbol{\sigma}d\Omega = \int_{\partial\Omega_N} <\textbf{v}, \textbf{f}_s>dA + \int_{\Omega} <\textbf{v}, \textbf{f}_{b}>d\Omega\\
\Rightarrow &-\int_{\Omega} <\boldsymbol{\phi}_j, \rho\sum_{i}\ddot{u}_i\boldsymbol{\phi}_i>d\Omega + \int_{\Omega}\nabla\boldsymbol{\phi}_j \colon \boldsymbol{\sigma}d\Omega = \int_{\partial\Omega_N} <\boldsymbol{\phi}_j, \textbf{f}_s>dA + \int_{\Omega} <\boldsymbol{\phi}_j, \textbf{f}_{b}>d\Omega\\
\Rightarrow &-\int_{\Omega} \rho\sum_{i}\ddot{u}_i<\boldsymbol{\phi}_j, \boldsymbol{\phi}_i>d\Omega + \int_{\Omega} \nabla\boldsymbol{\phi}_j \colon \mathcal{C} \colon \textbf{L}(\sum_{i} u_i \boldsymbol{\phi}_i)d\Omega = \int_{\partial\Omega_N} <\boldsymbol{\phi}_j, \textbf{f}_s>dA + \int_{\Omega} <\boldsymbol{\phi}_j, \textbf{f}_{b}>d\Omega\\
\Rightarrow &-\int_{\Omega} \rho\sum_{i}\ddot{u}_i<\boldsymbol{\phi}_j, \boldsymbol{\phi}_i>d\Omega + \int_{\Omega} \nabla\boldsymbol{\phi}_j \colon \mathcal{C} \colon \sum_{i} u_i \textbf{L}(\boldsymbol{\phi}_i)d\Omega = \int_{\partial\Omega_N} <\boldsymbol{\phi}_j, \textbf{f}_s>dA + \int_{\Omega} <\boldsymbol{\phi}_j, \textbf{f}_{b}>d\Omega\\
\Rightarrow &-\int_{\Omega} \rho\sum_{i}\ddot{u}_i<\boldsymbol{\phi}_j, \boldsymbol{\phi}_i>d\Omega + \int_{\Omega} \sum_{i} u_i \nabla\boldsymbol{\phi}_j \colon \mathcal{C} \colon \textbf{L}(\boldsymbol{\phi}_i)d\Omega = \int_{\partial\Omega_N} <\boldsymbol{\phi}_j, \textbf{f}_s>dA + \int_{\Omega} <\boldsymbol{\phi}_j, \textbf{f}_{b}>d\Omega\\
\Rightarrow &-\sum_{i}\ddot{u}_i\int_{\Omega} \rho<\boldsymbol{\phi}_j, \boldsymbol{\phi}_i>d\Omega + \sum_{i} u_i \int_{\Omega} \nabla\boldsymbol{\phi}_j \colon \mathcal{C} \colon \textbf{L}(\boldsymbol{\phi}_i)d\Omega = \int_{\partial\Omega_N} <\boldsymbol{\phi}_j, \textbf{f}_s>dA + \int_{\Omega} <\boldsymbol{\phi}_j, \textbf{f}_{b}>d\Omega
\end{split}
\end{equation}

The above equality holds for all $\boldsymbol{\phi}_j$. Thus we get a linear system that looks like following

\begin{equation}\nonumber
\textbf{M} \ddot{\textbf{u}} + \textbf{K}\textbf{u} = \textbf{f}
\end{equation}

where $\textbf{M}_{i,j} = -\int_{\Omega} \rho<\boldsymbol{\phi}_j, \boldsymbol{\phi}_i>d\Omega$ and $\textbf{K}_{i,j} = \int_{\Omega} \nabla\boldsymbol{\phi}_j \colon \mathcal{C} \colon \textbf{L}(\boldsymbol{\phi}_i)d\Omega$. \textbf{f} is just the total force of both surface traction and body force.

Such derivation can also be found in most of the FEM solid mechanics books but from a minimal virutal work point of view, which basically formuate potential energy first and then take variation.

Finally, we are going to talk a little bit about the relation between our derivation and those popular in computer graphics paper. We take ~\cite{Muller:2004:IVM} as an example. In that paper, they first have

\begin{equation}\nonumber
\textbf{u}(\textbf{x}) = \textbf{H}_e(\textbf{x})\hat{\textbf{u}}
\end{equation}

Here $\textbf{H}_e(\cdot)$ is just picewise linear nodal basis function whose support is elementwise. $\hat{\textbf{u}}$ is the nodal weight which are displacement in this sense. Then they have

\begin{equation}\nonumber
\boldsymbol{\epsilon} = \textbf{B}_e\hat{\textbf{u}}
\end{equation}

Here $\textbf{B}_e$ is a linear strain mapping which is Cauchy strain tensor in that case, corresponding to geometric linearity. Next, they have

\begin{equation}\nonumber
\boldsymbol{\sigma} =\textbf{E}\cdot\boldsymbol{\epsilon}
\end{equation}

where \textbf{E} a second order elasticity tensor for isotropic material which is forth order tensor in our case dealing with even anisotropic cases.

Then they construct their stiffness matrix as

\begin{equation}
\textbf{K}_e = V_e\textbf{B}_e^T\textbf{EB}_e
\end{equation}

which is similar to our case if you take our linear differential mapping as gradient operator and $\textbf{B}_e$ as gradient of basis functions, which really is. Then the volume integral gives you the volume multiplier. They also mentioned about stiffness warping, whose incorpration is easy to understand. However, we are not going to talk about it here as it's not embeded into FEM framework. One can take it as a strategy or preprocessing for FEM to deal with nonlinearity.

Finally, we wanna point out our mass matrix is not supposed to be diagonal even with elementwise support basis function. However, one can make use of matrix lumping techniques to sparsify or diagonalize it sacrificing accuracy a little bit but gained efficiency especially when dealing with KKT matrix in constrained dynamics.